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Abstract 

We build upon some new ideas in direct transcrip- 
tion methods developed within the Advanced Con- 
cepts Team to introduce two improvements to the 
Sims-Flanagan transcription for low-thrust trajec- 
tories. The obtained new algorithm is able to pro- 
duce an operational trajectory accounting for the 
real spacecraft dynamics and adapting the segment 
duration on-line improving the final trajectory op- 
timality. 

1 Introduction 

A direct optimisation method proposed by Sims 
and Flanagan [12] suggests that low-thrust trajec- 
tories can be modeled as a series of impulsive AV 
connected by conic arcs. The method is fast and 
robust and has been applied in previous works for 
preliminary mission design [15, 16]. However with 
its impulsive AV transcription, the Sims-Flanagan 
method can fail to accurately represent the actual 
dynamical model unless the number of impulses is 
increased and thus at the cost of slowing down the 
overall optimisation. 

In this paper, we study new transcription meth- 
ods to improve the accuracy of the Sims and Flana- 
gan model without increasing the dimension of the 
problem. The works extends some of the ideas pre- 
sented at the fifth international meeting on celestial 
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mechanics CELMEC V [5]. The first improvement 
is to replace the impulses with continuous thrust, 
where low-thrust arcs are numerically propagated. 
The magnitude and direction of the thrust are part 
of the optimisation variables and are assumed to be 
constant throughout a segment. Perturbations can 
be included, in the propagation, to further improve 
the fidelity of the model. The modification intro- 
duces a performance penalty due to the higher com- 
putational costs of the integration with respect to 
a simple Keplerian propagation between impulses. 
In order to tackle this issue we introduce the use of 
Taylor integration [G] methods in place of the com- 
monly employed Runge-Kutta-Fehlberg scheme, re- 
ducing the performance loss by almost one order of 
magnitude. 

A second improvement we introduce is to allow 
the time mesh to be optimised together with the 
trajectory. While this is a long unsolved issue in di- 
rect method for trajectory optimisation, we manage 
to obtain an efficient algorithm by introducing the 
Sundman transformation [13], in which the indepen- 
dent variable is changed from time to s, where the 
time variation of s is inversely proportional to the 
radial distance. By doing so, and adding only one 
constraint to the optimisation problem, segments 
are automatically distributed more densely near the 
central body (where speed is usually higher) along 
the optimal solution and thus on-line mesh adap- 
tation is obtained at the cost of an acceptable per- 
formance loss. We present a numerical example to 
compare results between the original and the new 
methods. The resulting tool has the further ad- 
vantage of being suitable for different phases of the 



mission design, from preliminary, where global op- 
timisation methods need a rather simple and low- 
dimensional transcription, to operational where dy- 
namics need to be accounted for in a precise manner 
and optimality is sought. 

a Matchpoint 

• Planet 

• Segment midpoint 
— ► Impulsive AV 




Figure 1: Impulsive AV transcription of a low- 
thrust trajectory, after Sims and Flanagan [12] 



2 The Original Sims- Flanagan 
Model 

In 1999, Sims and Flanagan proposed a di- 
rect method for optimising low-thrust trajectories, 
where later the software packages GALLOP [8] and 
MALTO [11] are developed based on this model. 
Figure 1 briefly illustrates such a trajectory model. 
The whole trajectory is divided into legs which be- 
gin and end with a planet. Low-thrust arcs on each 
leg are modeled as sequences of impulsive maneu- 
vers AV, connected by conic arcs. We denote the 
number of impulses (which is the same as the num- 
ber of segments) with N. The AV at each segment 
of equal duration should not exceed a maximum 
magnitude, AV max , where AY max is the velocity 
change accumulated by the spacecraft when it is 
operated at full thrust during that segment: 

AV max = (F max /m)(T f - T )/N (1) 

where F max is the maximum thrust of the low- 
thrust engine, m is the mass of the spacecraft, 



To and Tf is the initial and final time of a leg. 
The spacecraft mass is propagated using the rocket 
equation [14]: 

m i+ i = nii exp(-AVi/g I S p) (2) 

where the subscript i denotes the mass and AV 
on the i-th segment, go is the standard gravity 
(9.80665 m/s 2 ), and I sp is the specific impulse of 
the low-thrust engine. 

At each leg, trajectory is propagated (with 
a two-body model) forward and backward to 
a matchpoint (usually halfway through a leg), 
where the spacecraft state vector becomes S m / = 
{r x ,r y , r z ,v x ,v y , v z ,m} mf (and similarly for S m6 ), 
where r and v are respectively the position and 
velocity of the spacecraft and the subscripts rep- 
resents the Cartesian x, y, z components. The 
forward- and backward-propagated half-legs should 
meet at the matchpoint, or the mismatch in posi- 
tion, velocity, and mass: 

S m /-S m6 = {Ar x , Ar y , Ar z , Av x , Av v , Av z , Am} 

(3) 

should be less than a tolerance in order to have a 
feasible trajectory. 

The problem is transcripted into a nonlinear pro- 
gramming problem (NLP), where the objective is to 
maximize the final spacecraft mass subjected to the 
constraints on the maximum AV and the state mis- 
match, while the decision variables of the problem 
are listed below: 

• the departure epoch To 

• the departure velocity relative to the earth 

• for each leg and each segment, the magnitude 
of the impulse and direction 

• for each swingby, the incoming and outgoing 
velocities relative to the planet 

• for each swingby j, the swingby epoch Tj 

• the arrival epoch Tf 

For a rendezvous mission, the arrival velocity to 
the destination is not included in the set of vari- 
ables, as it is, by construction of the model, zero rel- 
ative to the planet. To solve the NLP, we use a soft- 
ware package called SNOPT [4] [3], which imple- 
ments sequential quadratic programming (SQP). 
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tion, the optimization process becomes slower as 
the optimization relies on a numerical integration 
scheme of a higher complexity with respect to a 
simpler ballistic arc solver. Efficiency is essential to 
keep the CPU time penalty at a minimum level. For 
this purpose, we report the comparison, in terms 
of CPU time and accuracy, among classical Runge- 
Kutta-Fhelberg methods and the Taylor integra- 
tion [6]. The tests have been done having in mind 
the typical algorithm call done during an optimiza- 
tion procedure that uses our approach. An improv- 
ment of one order of magnitude in CPU time (while 
keeping the integration accuracy to the same level) 
is found. 



Figure 2: Schematics of the three trajectory mod- 
els. Top: Impulisve AV; Middle: Continuous 
thrust is the time-space; Bottom: Continuous 
thrust in the s-space 

3 Improvement to Trajectory 
Model 

From the beginning, the use of the Sims-Flanagan 
model is limited to preliminary mission design, in 
which the results are not expected to be accurate 
up to the operation level. In terms of the fidelity of 
the model, there are two areas in the Sim-Flanagan 
model that can lead to loss in accuracy: (1) the 
use of impulses; and (2) the insufficient number of 
segments. To address these issues, we introduce 
two improvements: 

• Impulsive AV are replaced by continuous 
thrust to improve the fidelity of the trajec- 
tory dynamics. In order to keep a reasonably 
low computing time, we employ a Taylor inte- 
gration scheme showing an order of magnitude 
performance gain with respect to the classical 
Runge-Kutta methods. 

• An adaptive time mesh is obtained on the seg- 
ments via the Sundman transformation to im- 
prove on the optimality of the final trajectory. 

3.1 The Taylor Integration Method 

When replacing the original impulsive AV tran- 
scritpion with a continuous fixed thrust transcrip- 



Comparison set-up We consider the set of dif- 
ferential equations describing the motion of a space- 
craft subject to a fixed thrust force in the interplan- 
etary medium. This 'fixed thrust problem' is at 
the basis of the ideas on direct transcription meth- 
ods presented by some of these authors during the 
fifth international meeting on celestial mechanics 
CELMEC V ([5]) and that motivated the current 
paper. Since in the proposed new direct transcrip- 
tion method the fixed thrust problem needs to be 
solved a large amount of times (in each segment) 
and with diverse initial conditions and thrust vec- 
tors we focus the algorithmic comparison to those 
cases representative of such a process. The equa- 
tions, in a non dimensional form, are the following: 



r = \f'- 

X\ = X4 
X 2 = X 5 
X 3 = X 6 
&4 
X 5 
X 6 



xi/r 3 
x 2 /r 3 
x 3 /r 3 



(4) 
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The mass is not considered for the purpose of 
this comparison, but will be included in the tra- 
jectory model described later. To test the integra- 
tion schemes, N = 10000 different Cauchy prob- 
lems have been generated at random considering 
Xi(0),Ui uniformly distributed in Xj(0) G [0.1,2] 
and u t £ [0.0001,0.01]. The final integration 
time has been also set to be random and tf <G 
[7r/20, IOtt]. The same problems where solved using 
a Runge-Kutta-Fhelberg integration scheme (in the 
implementation of the GAL libraries [2] ) and a Tay- 
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lor integration scheme (implemented using the tool 
"taylor" [7]). In order to test the speed and the 
precision of the solvers, we propagate each prob- 
lem from Xi(0) for tf, we then take the result and 
propagate backwards for tf reaching the point x{ . 
By doing this, as we know the exact result of the 
propagation that is Xi(0), we evaluate the preci- 
sion of the propagation defining the propagation er- 
ror as err — Y^h=i ( x i(ty ~ x {j ■ Each algorithm 
is tested on the same set of randomly generated 
Cauchy problems. In all cases, no minimum step 
size is used and the same parameter e is passed to 
the RKF integrators as the absolute error, and to 
the Taylor integrator as both absolute and relative 
error. The initial trial stepsize of 0.1 is set to the 
RKF integrators. 



Results From the results oulined in Table 1 it is 
clear that the Taylor integrator is outperforming 
the RKF both in speed and accuracy confirming 
in the low-thrust fixed direction problem the same 
performance gain levels already reported in past 
literature [6, 10]. From the table we may also, em- 
pirically, establish that e = 10 -10 is a good compro- 
mise between speed and accuracy and can thus be 
used as a default parameter to call the Taylor inte- 
grator. The speed gained by employing the Taylor 
integrator is roughly one order of magnitude. 



3.2 The Sundman Transformation 

In his celebrated paper [13] Karl Sundman intro- 
duces a simple differential transformation for the 
time variable, to regularize the otherwise singular 
three body problem. The Sundman transformation 
dilates the time metric introducing a new variable 
s defined through the relation ds — dt/r guaran- 
teeing an asymptotically slower flow near the sin- 
gularities. In a trajectory propogation this same 
property turns out to be quite useful if equally 
spaced segments are considered in the s domain 
rather than in the t domain. Let us for example 
consider Eq.(4) and use the Sundman transforma- 
tion, we obtain the following set of equations: 



r = yjx\ + x\ + x\ 
x\ = x^r 
±2 = x 5 r 
x 3 = x 6 r 

±4 = —xi/r 2 + u\r 
&5 = -x 2 /r 2 + u 2 r 
%6 = -~x 3 /r 2 + u 3 r 
i = r 

To demonstrate the effect of such a transforma- 
tion on a numerical mesh, we take a circular orbit 
of radius one and propagate it forward with a con- 
stant thrust aligned along the x axis. In Figure 
3a we visualize the obtained orbit using a uniform 
sampling in time, while in Figure 4b the same tra- 
jectory is visulaized using the same number of sam- 
ples, but equally spaced in the Sundman variable 
domain. The same is done in Figures 4a-b for a 
constant tangential thrust. These pictures clearly 
show the problem with using points equally spaced 
in time to define a mesh for a numerical algorithm: 
due to the conservation of energy the closer we get 
to the singularity the more potential energy we lose 
and thus acquire in terms in kinetic energy. This 
pumps up the body velocity substantially creating 
an unequal distribution of segments length bound 
to create numercial difficulties. The use of the s 
variable is one of the possible transformations able 
to alleviate such a problem. 

4 New trajectory models 

The implementation of the ideas reported above 
leads to new trajectory models that, when opti- 
mized, result in significant improvments on the 
optimality and feasibility over the original Sims- 
Flanagan method. 

4.1 Continuous thrust time-space 
propagation 

In this model the impulses at each segment are re- 
placed with continouous thrusts (F x ,F y ,F z ) which 
are assumed to be constant within the segment (see 
the middle scheme in Figure 2). Each leg of the tra- 
jectory is propagated forward and backward with 
equal-duration segments as before. The propaga- 
tion of the trajectory changes from pure Keplerian 
to integration of the ordinary differential equations: 
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Table 1: Algorithm Performance. Speed is measured in seconds and refers to all the N integrations 
(forward+backward) . Max. Err. is the maximum integration error made in the N propagations. Note 
that this is a real error, not an estimation as explained above. 
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r z = v z 

v x = ~nr x /r 3 + Fx/m 
v y = -nr y /r 3 + Fy/m 
v z = —/ir z /r 3 + F z /m 
m = -F/(g I sp ) 



4.2 Continuous thrust s-space prop- 
agation 

For the continuous thrust s-space method, we ap- 
ply the Sundman transformation [13] to change the 
independent variable from time t to s,and the dif- 
ferential equations in s-space becomes: 



r - 


= ^rl + rl + 


^ z 
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= \j F t + F2 y 


+ F 2 

1 z 


fx 


= rv x 




*V 


= rv y 




fz 


= rv z 




V x 


= -nr x /r 2 -f 


rFx/m 


Vy 


= -^r y /r 2 + 
= -/ir z /r 2 + 


rFy/m 


V z 


rF z /m 


771 


= - rF l (go I s P ) 


i = 


- r 





where the derivatives are here meant to be taken 
with respect to the independent variable s. Here 
each leg of the trajectory is propagated forward 
from So and backward from s/ in equal s-space 
(As). The time between the mesh is no longer con- 
stant and it is proportional to the radial distance r, 
which implies a shorter time mesh (a finer grid size 
or segment) is used when the spacecraft is closer 
to the central body. The 8th differential equation 
gives the condition for matching the time difference 
between the two endpoints: 

T f -T = f ' ' rds (5) 

J s 
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Figure 3: A trajectory sampled with the same number of equally a) time spaced segments b) s-spaced 
segments 





b) 

Figure 4: A trajectory sampled with the same number of equally a) time spaced segments b) s-spaced 
segments 



To implement the s -space method for optimisa- 
tion, we assume So to be zero and solve for st to 
satisfy Eq. 5, which means an additional variable 
and constraint are added for each leg. 



5 Numerical Example 

We demonstrate our new methods with a sample 
low-thrust mission to Mercury. Table 2a summa- 
rizes the mission specification for a spacecraft sim- 
ilar to the Deep Space 1 mission [9]. Instead of 
using a model of the SEP (solar electric propul- 
sion) system, we simplify the problem a bit here by 



assuming a constant thrust and constant specific 
impulse engine. The launch and arrival dates are 
kept frozen for the test. 

Figure 5a shows the trajectory found by the im- 
pulsive model, where the black dots denote the 
midpoint of the segments and the red lines repre- 
sent the impulses. It is visually clear that for this 
fast rotating trajectory, the impulsive propagation 
method might not be able to have a fair represen- 
tation of the actual low-thrust trajectory. With 
the same number of segments (30), the continuous 
thrust time-space model is able to fill up the " gap" 
between the impulses with low-thrust arcs (shown 
as red curves in figure 5b). Even with a lower final 
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Parameters 



Values 



Initial mass of the spacecraft 

Maximum thrust 

Specific impulse 

Launch date 

Arrival date 

Launch I^o 

Time of flight 



660.0 kg 
92.3 mN 
3,337 s 
Apr. 9, 2007 
Aug. 22, 2013 
< 2.0 km/s 
6.37 years 



Method 




No. of 


Optimal Final 






Variables 


Mass, kg 


Impulsive 




97 


392.9 


Continuous t- 


space 


97 


382.4 


Continuous s 


-space 


98 


387.0 



(b) Optimal Final Mass for the Earth-Mercury Mission. Note 
that only the continuous trajecctories are feasible and thus can 
be used up to late design phases. d 



(a) Parameters for an Earth-Mercury Rendezvous Mission. 

Table 2: Trajectories parameters and results. 




mass (see Table 2b), the trajectory found by the 
continuous thrust method satisfies the real dynam- 
ical model which can be used as a trajectory for 
the actual mission (after adding other perturbative 
forces). The trajectory modelling can be further 
improved by converting to the s-space propagation 
method in figure 6a. In the s-space method, we 
note that when the spacecraft is near Mercury, the 
grid size is smaller than it is near Earth's orbit, 
while in the time-space method the segment dura- 
tion is constant (see figure 6b). In this example, the 
implementation of the s-space propagation can au- 
tomatically adapt the grid size (segment duration) 
to fit the radial distance (and hence the speed of 
the spacecraft) of the trajectory during the optimi- 
sation. A smarter choice of the segment duration 



along the trajectory allows the spacecraft to up- 
date its control (i.e. thrust) and therefore be more 
efficient, which explains why the final mass of the 
s-space is higher than the time-space method. 

6 Conclusions and Future 
Work 

We have successfully extended the Sims-Flanagan 
model to include the full dynamics of low-thrust 
trajectory. The change of independent variable via 
Sundman transformation can further improve the 
results through online adaptive time-mesh during 
the optimisation. In the future, we hope to inves- 
tigate a more general form of the Sundman trans- 
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Figure 6: Trajectory found using the Sundman transformation. 



formation [1] and to perform some benchmarking 
of the new methods to compare their convergence 
speed and the accuracy. 
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